#include "types.h"
#include "state_vector.h"
#include "convective.h"

//------------------------------------------------------------------------------

// calc r_i

double calc_r(double i) {

  return i * params->dr;

}

//------------------------------------------------------------------------------

double calc_vr_center_average(sv_t *y, int i, int j, int k) {

  return 0.5 * ( sv_get_vr(y, i, j, k) + sv_get_vr(y, i+1, j, k) );

}

//------------------------------------------------------------------------------

double calc_vt_center_average(sv_t *y, int i, int j, int k) {

  return 0.5 * ( sv_get_vt(y, i, j, k) + sv_get_vt(y, i, j+1, k) );

}

//------------------------------------------------------------------------------

double calc_vz_center_average(sv_t *y, int i, int j, int k) {

  return 0.5 * ( sv_get_vz(y, i, j, k) + sv_get_vz(y, i, j, k+1) );

}

//------------------------------------------------------------------------------

double calc_vr_edge_t_average(sv_t *y, int i, int j, int k) {

  return 0.5 * ( sv_get_vr(y, i, j, k) + sv_get_vr(y, i, j, k-1) );

}

//------------------------------------------------------------------------------

double calc_vz_edge_t_average(sv_t *y, int i, int j, int k) {

  return 0.5 * ( sv_get_vz(y, i, j, k) + sv_get_vz(y, i-1, j, k) );

}

//------------------------------------------------------------------------------

double calc_vr_edge_z_average(sv_t *y, int i, int j, int k) {

  return 0.5 * ( sv_get_vr(y, i, j, k) + sv_get_vr(y, i, j-1, k) );

}

//------------------------------------------------------------------------------

double calc_vt_edge_z_average(sv_t *y, int i, int j, int k) {

  return 0.5 * ( sv_get_vt(y, i, j, k) + sv_get_vt(y, i-1, j, k) );

}

//------------------------------------------------------------------------------

double calc_vz_edge_r_average(sv_t *y, int i, int j, int k) {

  return 0.5 * ( sv_get_vz(y, i, j, k) + sv_get_vz(y, i, j-1, k) );

}

//------------------------------------------------------------------------------

double calc_vt_edge_r_average(sv_t *y, int i, int j, int k) {

  return 0.5 * ( sv_get_vt(y, i, j, k) + sv_get_vt(y, i, j, k-1) );

}

//------------------------------------------------------------------------------

double convective_central_r(sv_t *y, int i, int j, int k) {

  double temp;

  temp = - 1.0 / calc_r(i) / params->dr * ( calc_r(i+0.5) * calc_vr_center_average(y, i, j, k) * calc_vr_center_average(y, i, j, k) - calc_r(i-0.5) * calc_vr_center_average(y, i-1, j, k) * calc_vr_center_average(y, i-1, j, k) );

  temp -= 1.0 / calc_r(i) / params->dtheta * ( calc_vt_edge_z_average(y, i, j+1, k) * calc_vr_edge_z_average(y, i, j+1, k) - calc_vt_edge_z_average(y, i, j, k) * calc_vr_edge_z_average(y, i, j, k) );

  temp -= 1.0 / params->dz * ( calc_vz_edge_t_average(y, i, j, k+1) * calc_vr_edge_t_average(y, i, j, k+1) - calc_vz_edge_t_average(y, i, j, k) * calc_vr_edge_t_average(y, i, j, k) );

  return temp;

}

//------------------------------------------------------------------------------

double convective_central_t(sv_t *y, int i, int j, int k) {

  double temp;

  temp = - 1.0 / calc_r(i+0.5) / params->dr * ( calc_r(i+1.0) * calc_vr_edge_z_average(y, i+1, j, k) * calc_vt_edge_z_average(y, i+1, j, k) - calc_r(i) * calc_vr_edge_z_average(y, i, j, k) * calc_vt_edge_z_average(y, i, j, k) );

  temp -= 1.0 / calc_r(i+0.5) / params->dtheta * ( calc_vt_center_average(y, i, j, k) * calc_vt_center_average(y, i, j, k) - calc_vt_center_average(y, i, j-1, k) * calc_vt_center_average(y, i, j-1, k) );

  temp -= 1.0 / params->dz * ( calc_vz_edge_r_average(y, i, j, k+1) * calc_vt_edge_r_average(y, i, j, k+1) - calc_vz_edge_r_average(y, i, j, k) * calc_vt_edge_r_average(y, i, j, k) );

  return temp;

}

//------------------------------------------------------------------------------

double convective_central_z(sv_t *y, int i, int j, int k) {

  double temp;

  temp = - 1.0 / calc_r(i+0.5) / params->dr * ( calc_r(i+1.0) * calc_vr_edge_t_average(y, i+1, j, k) * calc_vz_edge_t_average(y, i+1, j, k) - calc_r(i) * calc_vr_edge_t_average(y, i, j, k) * calc_vz_edge_t_average(y, i, j, k) );

  temp -= 1.0 / calc_r(i+0.5) / params->dtheta * ( calc_vt_edge_r_average(y, i, j+1, k) * calc_vz_edge_r_average(y, i, j+1, k) - calc_vt_edge_r_average(y, i, j, k) * calc_vz_edge_r_average(y, i, j, k) );

  temp -= 1.0 / params->dz * ( calc_vz_center_average(y, i, j, k) * calc_vz_center_average(y, i, j, k) - calc_vz_center_average(y, i, j, k-1) * calc_vz_center_average(y, i, j, k-1) );

  return temp;

}

//------------------------------------------------------------------------------

